# Estimate our main model & define the weight matrix within Province.

# remove everything
rm(list=ls())

# =================================================================================
# ---------------------------------------------------------------------
# Drop outliers in dataset "Data" according to the "v"th column
# at each percentage "p" at bottom and top. The default value of
# p is 0.005. v is a vector of numbers.
# ---------------------------------------------------------------------

outlier <- function(Data, v, p=0.005) {
  out <- c()
  for(i in v) {
    out <- rbind(out, quantile(Data[,i], probs=c(p, 1-p)))
  }
  j <- 1
  for(i in v) {
    Data <- Data[Data[,i]>out[j,1] & Data[,i]<out[j,2], ]; print(dim(Data))
    j <- j+1
  }
  return(Data)
}

# =================================================================================

# Load packages
library(foreign)
library(gmm)
library(plm)

# read data industry-by-industry
Data <- read.dta('Ind39.dta')

Data <- outlier(Data, which(names(Data)=="sm"), p=0.005)

# Delete obs according to variable "Foreign"
Data <- subset(Data, F1>=0 & F1<=1 & !is.na(Province) & sm<=1)
summary(Data)

# =============================================
# Step One
h.lnbe <- mean(log(Data$sm)); h.eta <- - log(Data$sm) + h.lnbe; h.e <- mean(exp(h.eta))
h.betam <- exp(h.lnbe)/h.e

# =============================================
# Step Two
# Generate weights & variables related to weights
for(i in 1:nrow(Data)) {
  s1.t <- (Data$F1 == 0)*(Data$Year==Data$Year[i])*(Data$Province==Data$Province[i]); s1.t[i] <- 0
  s1.b <- sum(s1.t)
  s1 <- if(s1.b>0) s1.t/s1.b else rep(0, length(s1.t));
  s2.t <- (Data$F1 > 0)*(Data$Year==Data$Year[i])*(Data$Province==Data$Province[i]); s2.t[i] <- 0
  s2.b <- sum(s2.t)
  s2 <- if(s2.b>0) s2.t/s2.b else rep(0, length(s2.t))
  Data$w1k1[i] <- sum(s1*Data$k1); Data$w1m1[i] <- sum(s1*Data$m1, na.rm = TRUE); Data$w1l1[i] <- sum(s1*Data$l1, na.rm = TRUE)
  Data$w2k1[i] <- sum(s2*Data$k1); Data$w2m1[i] <- sum(s2*Data$m1, na.rm = TRUE); Data$w2l1[i] <- sum(s2*Data$l1, na.rm = TRUE)
}

# Generate year dummies
TD <- model.matrix(~as.factor(Data$Year)-1); TD1 <- model.matrix(~as.factor(Data$Year-1)-1)[,-1]

# Define functions: "wphi" and "phi"
w1phi <- function(beta) {
  (1-h.betam)*Data$w1m1 - h.lnbe - log(Data$ppi1/Data$ppii1) - beta[1]*Data$w1k1 - beta[2]*Data$w1l1 - TD1%*%beta[3:10]
}

w2phi <- function(beta) {
  (1-h.betam)*Data$w2m1 - h.lnbe - log(Data$ppi1/Data$ppii1) - beta[1]*Data$w2k1 - beta[2]*Data$w2l1 - TD1%*%beta[3:10]
}

phi <- function(beta) {
  (1-h.betam)*Data$m1 - h.lnbe - log(Data$ppi1/Data$ppii1) - beta[1]*Data$k1 - beta[2]*Data$l1 - TD1%*%beta[3:10]
}


# -----------------------------------------------
# Optimization

objfn <- function(beta){
  y.star <- Data$y - h.betam*Data$m - beta[1]*Data$k - beta[2]*Data$l - TD%*%beta[3:11]
  p1 <- as.vector(phi(beta)); p2 <- as.vector(w1phi(beta)); p3 <- as.vector(w2phi(beta))
  var.poly <- poly(cbind(p1, Data$F1, p2, p3), degree = 2, raw = TRUE)
  mod1 <- lm(as.vector(y.star)~var.poly)
  OBJ <- sum(mod1$residuals^2)
  return(OBJ)
}

set.seed(1)
initial <- c(0.05, 0.1, 0.05, 0,0,0,0,0,0,0,0,0)
target <- optim(initial, objfn, method = "Nelder-Mead"); target

# results - elas
h.betak <- target$par[1]; h.betal <- target$par[2]

G.coef <- function(beta){
  y.star <- Data$y - h.betam*Data$m - beta[1]*Data$k - beta[2]*Data$l - TD%*%beta[3:11]
  p1 <- as.vector(phi(beta)); p2 <- as.vector(w1phi(beta)); p3 <- as.vector(w2phi(beta)) 
  var.poly <- poly(cbind(p1, Data$F1, p2, p3), degree = 2, raw = TRUE)
  mod1 <- lm(as.vector(y.star)~var.poly)
  return(mod1)
}
# results - partial effects
coef <- G.coef(target$par)$coef[2:15]
w1 <- as.vector(phi(target$par)); s1w1 <- as.vector(w1phi(target$par)); s2w1 <- as.vector(w2phi(target$par))

p.w  <- coef[1] + 2*coef[2]*w1 + coef[4]*Data$F1 + coef[7]*s1w1 + coef[11]*s2w1
p.F  <- coef[3] + coef[4]*w1 + 2*coef[5]*Data$F1 + coef[8]*s1w1 + coef[12]*s2w1
p.s1w1 <- coef[6] + coef[7]*w1 + coef[8]*Data$F1 + 2*coef[9]*s1w1 + coef[13]*s2w1
p.s2w1 <- coef[10] + coef[11]*w1 + coef[12]*Data$F1 + coef[13]*s1w1 + 2*coef[14]*s2w1

# The estimated w and w1
Data$w <- Data$y - h.betak*Data$k - h.betal*Data$l - h.betam*Data$m - TD%*%target$par[3:11] - h.eta
Data$w1 <- w1; Data$s1w1 <- s1w1; Data$s2w1 <- s2w1; Data$eta <- h.eta

# TIL
D <- data.frame(Data[,c('Firm','Year','Province')], pw=p.w, pF=p.F, ps1w=p.s1w1, ps2w=p.s2w1)
D <- pdata.frame(D, index = c("Firm", "Year"))
D$pF1 <- lag(D$pF, 1); 
D$DF1 <- Data$F1 == 0

D$temp1 <- ave(D$pF1*D$DF1, D$Year, D$Province, FUN = function(x) sum(x, na.rm = TRUE))
D$temp1[!is.na(D$pF1)] <- D$temp1[!is.na(D$pF1)] - (D$pF1*D$DF1)[!is.na(D$pF1)]
D$temp2 <- ave((!is.na(D$pF1))*D$DF1, D$Year, D$Province, FUN = function(x) sum(x, na.rm = TRUE))
D$temp2[!is.na(D$pF1)] <- D$temp2[!is.na(D$pF1)] - 1
D$sp=D$temp1/D$temp2
p.til1 <- as.numeric(D$ps1w*D$sp); 

D$temp1 <- ave(D$pF1*(!D$DF1), D$Year, D$Province, FUN = function(x) sum(x, na.rm = TRUE))
D$temp1[!is.na(D$pF1)] <- D$temp1[!is.na(D$pF1)] - (D$pF1*(!D$DF1))[!is.na(D$pF1)]
D$temp2 <- ave((!is.na(D$pF1))*(!D$DF1), D$Year, D$Province, FUN = function(x) sum(x, na.rm = TRUE))
D$temp2[!is.na(D$pF1) & D$temp2>0] <- D$temp2[!is.na(D$pF1) & D$temp2>0] - 1
D$sp=D$temp1/D$temp2; D$sp[abs(D$sp)==Inf]<-NA
p.til2 <- as.numeric(D$ps2w*D$sp); 

summary(p.s1w1); summary(p.s2w1)

# ===================================================================================
# ===================================================================================
# Section Two -- Bootstrap

# Calculate residuals and fitted values
resid <- G.coef(target$par)$residuals # ---------------
fit <- G.coef(target$par)$fitted.values # ---------------

# Start the Bootstrap
set.seed(123)
R <- 399

Bpw <- p.w; BpF <- p.F; Bps1w <- p.s1w1; Bps2w <- p.s2w1; Bptil1 <- p.til1; Bptil2 <- p.til2
Bbeta <- c(h.betak, h.betal, h.betam)
Bcoef <- coef

for(B in 1:R) {
  # Wild Bootstrap
  unif <- runif(length(resid),0,1)
  a <- (1-sqrt(5))/2
  prob.a <- (sqrt(5)+1)/(2*sqrt(5))
  b <- (1+sqrt(5))/2
  # ----------------------------------------
  # step 1
  h.eta.b <- ifelse(unif <= prob.a, a*h.eta, b*h.eta)
  lsm.b <- h.lnbe - h.eta.b
  
  h.lnbe.b <- mean(lsm.b); h.eta.b <- - lsm.b + h.lnbe.b; h.e.b <- mean(exp(h.eta.b))
  h.betam.b <- exp(h.lnbe.b)/h.e.b
  
  resid.b <- ifelse(unif <= prob.a, a*resid, b*resid)
  y.b <- h.betak*Data$k + h.betal*Data$l + h.betam.b*Data$m + fit + TD%*%target$par[3:11] + resid.b
  
  # ----------------------------------------
  # step 2
  
  objfn <- function(beta){
    y.star <- y.b - h.betam*Data$m - beta[1]*Data$k - beta[2]*Data$l - TD%*%beta[3:11]
    p1 <- as.vector(phi(beta)); p2 <- as.vector(w1phi(beta)); p3 <- as.vector(w2phi(beta))
    var.poly <- poly(cbind(p1, Data$F1, p2, p3), degree = 2, raw = TRUE)
    mod1 <- lm(as.vector(y.star)~var.poly)
    OBJ <- sum(mod1$residuals^2)
    return(OBJ)
  }
  target.b <- optim(initial, objfn, method = "Nelder-Mead")
  
  # results - elas
  h.betak.b <- target.b$par[1]; h.betal.b <- target.b$par[2]
  
  
  G.coef.b <- function(beta){
    y.star <- y.b - h.betam*Data$m - beta[1]*Data$k - beta[2]*Data$l - TD%*%beta[3:11]
    p1 <- as.vector(phi(beta)); p2 <- as.vector(w1phi(beta)); p3 <- as.vector(w2phi(beta))
    var.poly <- poly(cbind(p1, Data$F1, p2, p3), degree = 2, raw = TRUE)
    mod1 <- lm(as.vector(y.star)~var.poly)
    return(mod1)
  }
  # results - partial effects
  coef.b <- G.coef.b(target.b$par)$coef[2:15]
  w1.b <- as.vector(phi(target.b$par)); s1w1.b <- as.vector(w1phi(target.b$par)); s2w1.b <- as.vector(w2phi(target.b$par))
  
  p.w.b  <- coef.b[1] + 2*coef.b[2]*w1.b + coef.b[4]*Data$F1 + coef.b[7]*s1w1.b + coef.b[11]*s2w1.b
  p.F.b  <- coef.b[3] + coef.b[4]*w1.b + 2*coef.b[5]*Data$F1 + coef.b[8]*s1w1.b + coef.b[12]*s2w1.b
  p.s1w1.b <- coef.b[6] + coef.b[7]*w1.b + coef.b[8]*Data$F1 + 2*coef.b[9]*s1w1.b + coef.b[13]*s2w1.b
  p.s2w1.b <- coef.b[10] + coef.b[11]*w1.b + coef.b[12]*Data$F1 + coef.b[13]*s1w1.b + 2*coef.b[14]*s2w1.b
  
  # TIL
  D.b <- data.frame(Data[,c('Firm','Year','Province')], pw=p.w.b, pF=p.F.b, ps1w=p.s1w1.b, ps2w=p.s2w1.b)
  D <- pdata.frame(D.b, index = c("Firm", "Year"))
  D$pF1 <- lag(D$pF, 1); 
  D$DF1 <- Data$F1 == 0
  
  D$temp1 <- ave(D$pF1*D$DF1, D$Year, D$Province, FUN = function(x) sum(x, na.rm = TRUE))
  D$temp1[!is.na(D$pF1)] <- D$temp1[!is.na(D$pF1)] - (D$pF1*D$DF1)[!is.na(D$pF1)]
  D$temp2 <- ave((!is.na(D$pF1))*D$DF1, D$Year, D$Province, FUN = function(x) sum(x, na.rm = TRUE))
  D$temp2[!is.na(D$pF1)] <- D$temp2[!is.na(D$pF1)] - 1
  D$sp=D$temp1/D$temp2
  p.til1.b <- as.numeric(D$ps1w*D$sp); 
  
  D$temp1 <- ave(D$pF1*(!D$DF1), D$Year, D$Province, FUN = function(x) sum(x, na.rm = TRUE))
  D$temp1[!is.na(D$pF1)] <- D$temp1[!is.na(D$pF1)] - (D$pF1*(!D$DF1))[!is.na(D$pF1)]
  D$temp2 <- ave((!is.na(D$pF1))*(!D$DF1), D$Year, D$Province, FUN = function(x) sum(x, na.rm = TRUE))
  D$temp2[!is.na(D$pF1) & D$temp2>0] <- D$temp2[!is.na(D$pF1) & D$temp2>0] - 1
  D$sp=D$temp1/D$temp2; D$sp[abs(D$sp)==Inf]<-NA
  p.til2.b <- as.numeric(D$ps2w*D$sp); 
  
  # Store bootstrap results
  Bpw <- cbind(Bpw, p.w.b); BpF <- cbind(BpF, p.F.b); Bps1w <- cbind(Bps1w, p.s1w1.b); Bps2w <- cbind(Bps2w, p.s2w1.b);
  Bptil1 <- cbind(Bptil1, p.til1.b); Bptil2 <- cbind(Bptil2, p.til2.b)
  Bbeta <-cbind(Bbeta, c(h.betak.b, h.betal.b, h.betam.b))
  Bcoef <- cbind(Bcoef, coef.b)
  
}

apply(Bbeta, MARGIN = 1, sd)
summary(Bpw[,1]); apply(apply(Bpw, MARGIN = 2, summary), MARGIN = 1, sd)
summary(BpF[,1]); apply(apply(BpF, MARGIN = 2, summary), MARGIN = 1, sd)
summary(Bps1w[,1]); apply(apply(Bps1w, MARGIN = 2, summary), MARGIN = 1, sd)
summary(Bps2w[,1]); apply(apply(Bps2w, MARGIN = 2, summary), MARGIN = 1, sd)
summary(Bptil1[,1]); apply(apply(Bptil1, MARGIN = 2, summary), MARGIN = 1, sd)
summary(Bptil2[,1]); apply(apply(Bptil2, MARGIN = 2, summary), MARGIN = 1, sd)

save(BpF, Bpw, Bps1w, Bps2w, Bptil1, Bptil2, Bbeta, Bcoef, file = "Boot_ss_Dt.RData")

